#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <time.h>
#include <math.h>

using namespace std;

int N;
double Lx, Ly, Lz;
double dt, epsilon = 1, sigma = 1;
double K, P, E;
double max_v;

double *x, *y, *z;
double *x_prev, *y_prev, *z_prev;
double *vx, *vy, *vz;
double *ax, *ay, *az;

bool period_z_min = false;
bool period_z_max = false;

int iterations;
 
void initial_fromFile()
{
	cout << "dt: ";
	cin >> dt;

	cout << "the number of iterations: ";
	cin >> iterations;

	cout << "Lx: ";
	cin >> Lx;

	cout << "Ly: ";
	cin >> Ly;

	cout << "Lz: ";
	cin >> Lz;

	ifstream flIn;
	flIn.open("input.txt");

	flIn >> N;

	x = new double[N];
	y = new double[N];
	z = new double[N];

	x_prev = new double[N];
	y_prev = new double[N];
	z_prev = new double[N];

	vx = new double[N];
	vy = new double[N]; 
	vz = new double[N];

	ax = new double[N];
	ay = new double[N];
	az = new double[N];

	for( int i = 0; i<N; ++i )
	{
		flIn >> x[i];
		flIn >> y[i];
		flIn >> z[i];
		flIn >> vx[i];
		flIn >> vy[i];
		flIn >> vz[i];

		ax[i] = 0;
		ay[i] = 0;
		az[i] = 0;
	}

	flIn.close();
}

void destruction()
{
	delete[] x;
	delete[] y;
	delete[] z;

	delete[] x_prev;
	delete[] y_prev;
	delete[] z_prev;

	delete[] vx;
	delete[] vy;
	delete[] vz;

	delete[] ax;
	delete[] ay;
	delete[] az;
}

void periodic( double &x, double &y, double &z, double &vz)
{
	int lx;
	double mx;
	mx = frexp(Lx, &lx);

	int ly;
	double my;
	my = frexp(Ly, &ly);

	while( x < 0 || x > Lx || y < 0 || y > Ly || z < 0 || z > Lz)
	{
		int kx;
		double Mx;
		Mx = frexp(x, &kx);

		int ky;
		double My;
		My = frexp(y, &ky);

		int stxx = kx - lx;
		int styy = ky - ly;

		long double xx = pow((double)2,stxx);
		long double yy = pow((double)2,styy);

		if ( x < 0 )
		{
			if ( xx > 1 )
				x += Lx*xx;
			else
				x += Lx;
		}

		if ( x > Lx )
		{
			if ( xx > 1 )
				x -= Lx*xx;
			else
				x -= Lx;
		}

		if ( y < 0 )
		{
			if ( yy > 1 )
				y += Ly*yy;
			else
				y += Ly;
		}
		if ( y > Ly )
		{
			if ( yy > 1 )
				y -= Ly*yy;
			else
				y -= Ly;
		}

			if ( z<=0 )
		{
			z = -z; 
			vz = -vz;
			period_z_min = true;
		}

		if ( z>= Lz )
		{
			z = Lz - (z - Lz);
			vz = -vz;
			period_z_max = true;
		}
	}
}

void f( double r, double &force, double &potential )
{
	double ri = 1/r;
	double ri3 = ri*ri*ri;
	double ri6 = ri3*ri3;

	force = 48*epsilon*( pow(sigma,12)*ri6 - pow(sigma,6)/2 )*ri6*ri*ri;
	potential = 4*epsilon*ri6*( ri6* pow(sigma,12) - pow(sigma,6) );
}

void accel()
{
	double dx, dy, dz;
	double force, potential;
//	P = 0;

	for( int i = 0; i < N; i++ )
	{
		ax[i] = 0;
		ay[i] = 0;
		az[i] = 0;
	}

	for( int i = 0; i < N - 1; ++i )
	{
		for( int j = i + 1; j < N; ++j )
		{
			dx = x[i] - x[j];
			dy = y[i] - y[j];
			dz = z[i] - z[j];

			if ( fabs(dx) > .5*Lx )
			{
				if ( dx > 0 )
					dx = dx - Lx;
				else
					dx = dx + Lx;
			}

			if ( fabs(dy) > .5*Ly )
			{
				if ( dy > 0 )
					dy = dy - Ly;
				else
					dy = dy + Ly;
			}

			double r = sqrt( dx*dx + dy*dy + dz*dz );

			force = potential = 0;

			if( r == 0 )
				r = 0.01;

			f( r, force, potential );

			ax[i] = ax[i] + force*dx;
			ay[i] = ay[i] + force*dy;
			az[i] = az[i] + force*dz;

			ax[j] = ax[j] - force*dx;
			ay[j] = ay[j] - force*dy;
			az[j] = az[j] - force*dz;

//			P = P + potential;
		}
	}
}

void verlet()
{
//	K = 0;
	double *x_buf, *y_buf, *z_buf;
	double x_new, y_new, z_new;

	for ( int i = 0; i < N; i++ )
	{
		x_new = 2 * x[i] - x_prev[i] + ax[i] * dt * dt;
		y_new = 2 * y[i] - y_prev[i] + ay[i] * dt * dt;
		z_new = 2 * z[i] - z_prev[i] + az[i] * dt * dt;

		periodic(x_new, y_new, z_new, vz[i]);

		x_prev[i] = x_new;
		y_prev[i] = y_new;
		z_prev[i] = z_new;

		if (period_z_min == true)
		{
			z[i] = -z[i];
			period_z_min = false;
		}
		else
		if (period_z_max == true)
		{
			z[i] = Lz - (z[i] - Lz);
			period_z_max = false;
		}

	}

	x_buf = x_prev;
	x_prev = x;
	x = x_buf;

	y_buf = y_prev;
	y_prev = y;
	y = y_buf;

	z_buf = z_prev;
	z_prev = z;
	z = z_buf;

	for( int i = 0; i < N; ++i )
	{
		vx[i] = vx[i] + 0.5*ax[i]*dt;
		vy[i] = vy[i] + 0.5*ay[i]*dt;
		vz[i] = vz[i] + 0.5*az[i]*dt;
	}

	accel();

	for( int i = 0; i < N; ++i )
	{
		vx[i] = vx[i] + 0.5*ax[i]*dt;
		vy[i] = vy[i] + 0.5*ay[i]*dt;
		vz[i] = vz[i] + 0.5*az[i]*dt;


//		K = K + .5*( vx[i]*vx[i] + vy[i]*vy[i] );
	}
}

void eiler()
{
	double *x_buf, *y_buf, *z_buf;
	double x_new, y_new, z_new;

	x_buf = x;
	y_buf = y;
	z_buf = z;

	x = x_prev;
	y = y_prev;
	z = z_prev;

	x_prev = x_buf;
	y_prev = y_buf;
	z_prev = z_buf;

	for ( int i = 0; i < N; ++i )
	{
		x_new = x_prev[i] + vx[i]*dt;
		y_new = y_prev[i] + vy[i]*dt;
		z_new = z_prev[i] + vz[i]*dt;

		periodic(x_new, y_new, z_new, vz[i]);

		x[i] = x_new;
		y[i] = y_new;
		z[i] = z_new;

		if (period_z_min == true)
		{
			z_prev[i] = -z_prev[i];
			period_z_min = false;
		}
		else
		if (period_z_max == true)
		{
			z_prev[i] = Lz - (z_prev[i] - Lz);
			period_z_max = false;
		}

		vx[i] = vx[i] + 0.5*ax[i]*dt;
		vy[i] = vy[i] + 0.5*ay[i]*dt;
		vz[i] = vz[i] + 0.5*az[i]*dt;

	}

	accel();

	for( int i = 0; i < N; i++ )
	{
		vx[i] = vx[i] + 0.5*ax[i]*dt;
		vy[i] = vy[i] + 0.5*ay[i]*dt;
		vz[i] = vz[i] + 0.5*az[i]*dt;

//		K = K + .5*( vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
	}
}

void print_position(ofstream &fl)
{
	for (int i = 0; i < N; i++)
	{
		fl << x[i] << " " << y[i]<< " " << z[i] << "\n";
	}

	fl << "\n";
}

void print_energy(ofstream &fl)
{
	E = K + P;
	fl << "\n";
	fl << K << " " << P << " " << E;
	fl << "\n";
}

int main()
{
	initial_fromFile();

	/*P = 0;
	K = 0;
	E = 0;*/

//	print_position(fl);

	eiler();

//	print_position(fl);
//	print_energy(fl);

	for ( int i = 0; i < iterations; i++ )
	{
		verlet();
//		print_position(fl);
//		print_energy(fl);
	}

		ofstream flOut;
		flOut.open("output.txt");

		print_position(flOut);

		flOut.close();

	destruction();
	return 0;
}

